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DETERMINATION OF ORBIT OF A SPACECRAFT WITH RESPECT TO 
AN OBJECT IN A KNOWN CIRCULAR ORBIT 

By Robert Scott Dunning 
Langley Research Center 

SUMMARY 

A set of second-order differential equations of motion for a body in a planetocentric 
orbit has been derived and solved with the use of a cylindrical '’shell" coordinate system 
which describes the motion of an orbital vehicle as observed from another body in a 
known circular orbit about the planet. The solutions have been investigated for the spe- 
cific case of a body in orbit about the moon. The results of the study indicate improve- 
ment in accuracy over first-order theory with no serious adverse sensitivity to initial 
conditions and with applicability over a wide range of flight conditions. Orbit prediction 
is accurate over relatively long trajectories with large spatial separation of the vehicles. 
Some parts of the study are applicable to a manual apogee and perigee prediction scheme. 

INTRODUCTION 

In this paper, solutions are presented to a set of equations which describe the 
motion of a space vehicle in orbit about a planetary body as this motion would appear if 
seen from another space vehicle in orbit about the same planetary body. Thus, this 
report is concerned primarily with the analysis of the motion of a body in a moving rela- 
tive coordinate system. 

The method differs from the second -order method, which has been used success- 
fully in references 1 and 2, in that the coordinate system is chosen in such a manner as 
to eliminate or render negligible coupling between the various directions of translational 
motion in the mathematical treatment. The coordinate system used is similar to the shell 
system of reference 3 where, however, only first-order terms were considered. In this 
study, a cylindrical shell system has been adopted in preference to the spherical shell 
system of the latter reference. The method of attack is to examine the equations of 
motion of one vehicle as viewed from the frame of reference of the other vehicle by means 
of a set of second -order differential equations of motion. By making use of an exact first 
integral of one of the equations, solutions can be obtained that are accurate to second 
order. A simple perigee prediction method is also suggested which may be suitable for 
manual solution. 



A specific application of this analysis is the problem of trajectory prediction for a 
vehicle which is launched from a lunar orbit in an attempt to land on the lunar surface . 
The numerical computations are based on a circular orbit 200 kilometers above the lunar 
surface. The mathematical development is, however, thought to be quite general and, 
hence, is thought to be applicable to orbits about any spherical, gravitational body. 

SYMBOLS 

A integration constant for equation (17) defined in equation (20) 

a constant of integration which describes amplitude of departure from reference 

circular orbit 

c constant of integration determined by initial conditions, defined in equation (20) 

B,C,D,E, F constants determined by initial conditions defined in equation (20) 

H total energy of landing vehicle 

i,j,k orthogonal coordinates 

* 

x 

K constant, 1 

u>r s 

L Lagrangian, potential energy subtracted from kinetic energy 

M translation substitution used to solve equation (A30) and defined in 

equation (A36) 

m mass of landing vehicle 

r distance from center of planet to orbital vehicle 

r m radius of attracting body 

r distance from center of attracting body to reference orbit 

s 

t time 


2 



u 


dummy variable, a cos e 


V substitution variable used to solve equation (A30) 

X,Y,Z nondimensional coordinates of shell coordinate system centered on a body 
moving in a circular orbit about an attracting planet or satellite 

x, y,z dimensional coordinates of shell coordinate system centered on a body moving 

in a circular orbit about an attracting planet or satellite (see fig. 1) 

a,i 3, A functions of initial velocity component in x-direction defined in equation (20) 

A function of a power expansion in p defined in equation (20) 

e arbitrary constant of integration which describes epoch angle with respect 

to reference circular orbit 

6 angular coordinate in a cylindrical system with origin at center of attracting 

body and lying in plane of reference circular orbit, measured from positive 
x-direction clockwise 

y, T), k constant functions of initial out-of -plane velocity defined in equation (20) 

p gravitational constant 

p function of a, j3, and A defined in equation (20), hence a function of initial 

velocity in x-direction only 

r nondimensionalized or scaled time, cot 

<p ejection angle (see fig. 1) 

\f/ function of a power expansion in p defined in equation (20) 

co angular velocity 

directional angular velocity 
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Subscripts: 


max maximum 

o initial conditions 

p perigee 

Dots over symbols denote differentiation with respect to time or scaled time 

depending upon the section of the report in which they appear. Primes refer to first- 
order solution quantities. For constants, any consistent set of units may be used. In 
this paper the following values, which apply to the moon, were chosen: 

r m = 1736.5 kilometers 

r = r + 200 = 1936.5 kilometers 
s m 

/i= 4.8936 x 10^ meter^/second^ 
c v = 0.00082086 radians/second 
COORDINATES AND EQUATIONS OF MOTION 

Two assumptions are made about the physical nature of the problem: 

(1) The attracting planetary mass is a gravitational sphere. 

(2) The body upon which the coordinate system is centered is in a circular orbit. 

The coordinate system employed in this development is shown in figure 1 where 
r is the projection of a line connecting the center of the planet to the orbital vehicle upon 
the plane of the reference vehicle, z is a normal to this plane passing through the 
orbital vehicle, y is measured along r from the reference vehicle altitude to the pro- 
jection of the maneuvering vehicle with the positive direction upward, and x is measured 
in a curved arc backward along the flight path of the reference vehicle in the plane of the 
reference vehicle orbit to r. The coordinate system rotates about the origin with angu- 
lar velocity o>. 
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Direction of 
orbital motion 



Figure 1.- Coordinates employed in describing the motions of the vehicles. 
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This motion may be determined from the equation for the total energy and this 
method leads to an exact solution. However, this approach yields time as a function 
of y transcendentally and thus is not thought to be desirable since time normally is the 
available independent variable . 

The derivation of the equations of relative motion is given in detail in appendix A; 
therefore, only the more pertinent results are reviewed in this section. 

The equations of motion can be derived through the use of the Lagrangian which, for 
cylindrical coordinates, is 


L = imfr 2 + r 2 (0 - w) 2 + z 2 | + m — ^ (1) 

2 L J 6 * 7 ? 

Equation (1) is converted to shell coordinates (see appendix B) by means of the following 
substitutions: 


y + r g = r 

x = V 

z = z 


and use is made of the exact orbital expression for the angular velocity 


<o 2 = -l* 

r S 3 


to express the equations in the cylindrical shell system. Since the x-coordinate is cyclic, 
a first integral to its differential equation of motion is found immediately. The resulting 
differential equations of motion are 
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By solving equation (2) for 1 , 

<jOY 

x s 

it becomes possible to substitute the 
solution into equation (3) and hence 
remove all terms in x from the 
latter. Then, equation (3) is coupled 
to equation (4) only through the 
z-term. It is, however, found from 
experience that in orbits which are 
economically feasible, z is much 
smaller than y. Hence, there will 
be only a small amount of error in 


equation (3) if the term f — 

r 


is 


I78.4x)0 3 



assumed to be negligible even 
though this term is, strictly 
speaking, a second -order term in 
the mathematical sense. A spe- 
cific example of this relationship 
can be seen in figure 2 which shows 
the ratio of maximum z amplitude 
to maximum y amplitude for a 
typical lunar synchronous orbit over 
a range of economically feasible 
initial z velocity values. It is also 
found expedient to expand terms of 
the form when they occur. 

v\ n 

1 


I784xl0 3 



yp,max 


Figure 2.- Effect of maximum i excursions on perigee altitude tor 
synchronous orbits using exact equations. 


If the assumptions indicated are made, the equations which are to be solved are 


i_ + 2K£_ (K _ 1)= 3I^ 


cor c 


( 5 ) 


u 2 r s r s 


^ = _ ^ 


( 6 ) 


z | _z_ = 3yz 
w2r s r s r s 2 


(?) 
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where x Q , y Q , and z Q are arbitrary initial velocities and where 



1 


o? = 3K 2 - 2 
0 2 = K 2 - 1 
X = -6K 2 + 3 


In order to get an intuitive feel for the values of these constants, it is pointed out that 

x Q 

u>r s is the speed of the orbiting vehicle and hence that is never likely in a prac- 

wr s 

tical case to approach unity. In fact, for purposes of this report, it is assumed that x Q 

x 0 

is around 150 m/sec or less for the moon. Under these circumstances will be 

cor s 

less than 0.1. Thus, for a rough approximation as to order of magnitude, K is on the 
order of -1, c? is approximately unity, (ft is very small but can be either positive 
or negative, and X will be on the order of -3. 

The initial conditions assigned are those at the time of ejection when t = 0, 
x = y = z = 0, and x = x Q , y = y , z = z Q for all the cases studied in this report. 
Equation (6) is cast in the form of a one -dimensional anharmonic oscillator and hence 
can be solved with y as a time -dependent variable (ref. 4 and eq. (A31)). The solution 
of the equation of a one -dimensional anharmonic oscillator however is an elliptic integral 
and since this solution would prove difficult to incorporate in the solution of equations (5) 
and (7), a perturbation solution is accepted. This solution is then used in a similar man- 
ner to obtain an integral of equation (5) and, along with a first-order solution for equa- 
tion (7), can also be used to provide a second-order perturbation solution for equation (7). 
Hence, equations (5), (6), and (7) can be solved. A complete treatment of the solution of 
these equations is given in appendix A. 

The following solutions are to first order, where primes have been used to distin- 
guish the integration constants from those which will occur in second-order theory 



2 -| ■)cut - 2 Ka sin(cKi>t + e ') + sin € * 

a 2 ) « a 


( 8 ) 
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( 9 ) 


= a' cos(awt + e ') + 
r s or 


— = — — s i n w t 
r s wr s 


and where the integration constants a’ and e ' are given by 




€ ' = COS 


2 - 2 

, 4.^0 


Equations (11) and (12) must be used with a sign selection scheme. In order to do so, 
equation (9) can be expressed in a deterministic form without the epoch angle as 


v b 2 y Q B 2 

cos awt + sin owt + Sr 


Comparison of equations (13) and (9) gives 

B 2 

a' cos e ' = - Sr 

or 


a' sin e ' = - 


Hence, 


tan e * = 


9 



Since a is a frequency, which is a physical quantity, a negative a has no meaning. 
Hence, with a positive, the quadrant for e ' is determined by the respective signs in 
the numerator and denominator. Also /3 is strictly a function of x Q and always has 

2 x / x \ x 

the opposite sign, /3 = 1 2], since « 2. It then follows, once the quadrant 


a>r s \a> r s 


wr £ 


for 6 ' is known, that equations (14) and (15) give the sign for a' to be selected. It is 
found that a’ is always positive. These characteristics for all possible ejection quad- 
rants are given in table I. 

Solutions to the same differential equations but carried out to second order are 

x = r s |^A + Bo>t + C sin(po>t + e) + D sin 2(pu>t + e) + E sin 3(po>t + e) + F sin 4(pcot + e)j 

(17) 


y = r g ^ cos(pu>t + e) - ^ jl - | cos 2(pwt + e)J - 


(18) 


Z = IV 


tor 


— sin u>t + sin wt ^ cos(po>t + e) + cos 2(pwt + e) + - 
_4 4 4 


cos o>t sin (pwt + e ) + sin 2(po>t + e ) + 

|_2p 4p 2 2J 


(19) 


where the constant terms are given by 

*c 
’s 

,2 _ orr2 


K = — 1 
wr c 


a * = 3K* - 2 


/3 2 = K 2 - 1 


A = -6K a + 3 


p = (4A/3 2 + a 4 ) 1 ^ 4 = (-15K 4 + 24K 2 - 8) 


1/4 


(Equation continued on next page) 
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A = -(C sin e + D sin 2e + E sin 3e + F sin 4e) 


B = K 


c o (a 2 - p 2 ) ~ P 2 ) a 2 X 3o?a? 19 X 2 a 4 

T f ' * T T ““ T T ' . 


W r g K 


C = -K 


2a , 32c(c? - p 2 ) , 5Xa 

‘ • c 


4X 

3 


Xp 


2p c 


2p^ 


24 o 4 


T >-P'l a2 a 2 Q ! 2 Xa ! 2 X 2 a 4 

\P " 4p3 6p3 ” 4p3 j 


E = 


KXa c 

6p 3 


> ( 20 ) 


F = 


KX 2 a 4 

96p 5 


K = 


-3z 0 

Xa 2 

u>r s 

2p 2 

Xa 2 z. 


o 


Xa 2 (g 2 - p 2 ) 

a 


2X 


2cor s p^ 

3az„ 


y = 


o>r. 


A = 1 + p 2 + p 4 


.. p2 p4 

^ =1 + 4- + 16 


c = - sin e - 2^- sin 2e 
P 2p 


The determination of the integration constants, a and e , is somewhat more 
involved under second-order theory because of the necessity of solving a set of simul- 
taneous transcendental equations. Taking the y-equation and its first time derivative 
with the same end conditions as previously applied (namely, at t = 0, x = 0, y = 0, 
x = x Q , and y = y c ) and making the substitution u = a cos e yields the quartic equation 
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( 21 ) 



_2X1 U 4 + 3 + £ 


2 - «2\ 2 

*u + 


V 


2cP 


3p^ 


2X 


u + 


3 PV-4 



It is necessary to solve this equation numerically. Of the four possible values 
of u which satisfy this equation, only one will apply; this value is the value which 
approximates the expected amplitude of the y-expression. The computations made for 
this report give two imaginary roots and two real roots. The extraneous real root was 
found to be an order of magnitude larger than the correct root. Hence, comparison with 

first-order theory ^where it is found that u' = - affords a comparatively simple 

method of separating the roots. 

From equation (A58) and the definition of u, 


a = 


3py r 


-.2 


o>r 


,(3p 2 + 2Au) 


+ u 


( 22 ) 


and 


- „„ -1 u 

e = cos — 
a 


(23) 


It is assumed that approximately the same sign and quadrant selection scheme 
holds for the second-order integration constants as for the first-order integration con- 
stants, that is, a is always positive with the quadrant for e determined by <p. This 
relationship, however, is only approximately true. Examination of equation (21) for the 


special case where u = 0 gives for 


wr t 



(24) 


The terms on the right in equation (24) are functions of x Q only. Sketch (1) (see appen- 
dix A) shows how Xq is related to Y 0 when u = 0. This difference is the primary 
difference between first- and second-order theory. The sign on u along with the ejec- 
tion angle specifies the quadrant for e . This quadrant selection is summarized in 
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table II. Since the terms under the radical are all functions of 
angle is related to x Q by 


tan <p ~ - — = 




only, the ejection 


The parameters p 2 , c?, and X 2 are functions of x Q only; thus the quadrant 

of e changes depending upon cp, but no longer simply at ir/2, 7 r, 37r/2, and 2 tt as 

was true for the first-order solutions. Figure 3 is a plot of x Q as a function of <fi for 
u = 0 in dimensional terms which apply to the moon. The figure can be used therefore 
to determine for the known value of x Q and c p , the proper quadrant for e. As in first- 
order theory, a is always selected positive. 


aoxio 2 


i.q- 


*0. 

m/sec 


u = 0 


1 . 2 - 

in fourth quadrant^ 

.8 — 


• in third quadrant - 


40 


,\l 




- € in second quadrant 


u=0 


1 


in first quadrant*- 


80 120 


160 200 
<£, deg 


240 280 320 360 


Figure 3.- u=0 as a function of x Q and ejection angle 0. 

PHYSICAL INTERPRETATION OF TERMS 

The physical meaning of several of the terms in these solutions can be seen if the 
first-order solutions and second -order solutions are compared. 

y-Equation 

The most important from the standpoint of applications is the y-equation. It will be 
best to begin by comparing equation (9) (first-order solution) with equation (18) (second- 
order solution). In first-order theory, a sinusoidal term is added to a constant which is 
a function of initial velocity in the x-direction. Hence, the constant merely reflects a 
change in energy produced by ejection from the reference orbit and is a function of the 


13 



I I 


speeding up or slowing down of the orbital vehicle's angular rate. Superimposed upon 
this motion is an oscillatory motion above and below the mean altitude, the maximum 
amplitudes of which are apogee and perigee. In second-order theory, the same funda- 
mental characteristics appear. In this case, however, they are modified slightly to take 
better account of the inverse square nature of the planetary gravity field. The form is 
essentially the same as first-order theory presents, but the constant term is slightly 
modified to become the coefficient of the second term combined with the last term on the 
right in equation (18). In like manner, the oscillatory term is present although the coef- 
ficient is changed slightly. Hence, these two terms represent the same type of motion as 
was present in the first-order theory. Superimposed upon this term is an oscillatory 
term of twice the frequency of the primary term acting very much as a second harmonic. 
The primary effect of this term is to increase the apogee and to decrease the perigee. 
The function of this term is then to take account of the weakening of the restoring force 
as the gravity field decreases at longer distances from the center of the planet. 

Another way of viewing this phenomenon can be seen in figure 4 which shows the 
acceleration of the lander relative to the reference orbit in terms of exact first-order 
and second -order theory for a typical case. It can be seen that acceleration is a linear 
function of displacement over a reasonable altitude range under first-order theory. 



y. m 


Figure 4.- Acceleration of the orbital vehicle lander as a function of departure distance from the 
reference orbit for a typical case. 
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Under second-order theory, the acceleration is quadratic and forms a much better 
approximation to the exact conditions. In figure 4, /3 /a constitutes a change in mean 

equilibrium altitude and shows that, in general, the point of equilibrium is displaced, in 
this case downward. It is pointed out that this change will occur even for orbits of the 
same total energy as the reference orbit since orbits of the same period but different 
eccentricity will have a different mean altitude. Under second-order theory this mean 
displacement is defined a little more precisely by 


-r 


s 



+ 



(not shown in fig. 4). It can be seen that under the special condition that X = 0, this 
relation reduces to /3 ^/o? as in first-order theory. 


x -Equation 

Comparison of the x-equation (eq. (17)) with first-order theory (eq. (8)) shows that 
the first-order equation contains three main terms: a constant, a secular term, and a 
sinusoidal term. From the physics of the situation, it is clear that the phenomenon 
which was manifested as a constant mean separation in the y-equation (eq. (18)) is a 
time -dependent linear mean drift rate in the x-equation (eq. (17)). In general, the ejected 
vehicle would be expected to drift away over an interval of several orbital periods due to 
its difference in energy and angular momentum. For a synchronous orbit, of course, the 
coefficient of the secular term would be zero; thus, this term controls the rate of drift. 

In other respects the form of this equation is about the same as is obtained for the 
y-equation, that is, a series of higher order trigonometric terms. 


z -Equation 

In the z-equation, first-order theory predicts simple harmonic motion; so again 
there is a linear restoring force. This theory then assumes that altitude y does not 
have any appreciable effect upon the out-of -plane mode z. Although in trajectories of 
practical interest it is quite true that z does not have any marked effect upon y 
because the orbits are nearly coplanar, the converse of this statement is not true, and it 
is found that for orbits of moderately large eccentricity (for instance, on the order of 0.1), 
y-coupling into the z-equation can be very significant relative to the total, but admittedly 
small, amplitude in z. In fact, for the typical snychronous lunar orbit used for the 
numerical calculations in this paper, the coupled term in y and z amounts to one- 
third of the value of the pure z-term. Consequently, the second-order terms in equa- 
tion (19) represent coupling terms with y. 
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APOGEE AND PERIGEE PREDICTION 


The total orbital energy can be used in the following manner to predict apogee 
and perigee: 

Equation (5) is substituted into the orbital energy expression 

,2 



2 


H = i-m 
2 


.2 -2 

y + z + 


( y+r 4H 2 -7 =¥= 

' s ' I (y + r„) 2 + 


i/(y + r s) L + z ‘ 

The out-of -plane terms, z and z, are neglected since they are small. 
Use of equation (5) yields 

,2 


( 25 ) 


C = 


H 1 




t 2 . + 1 1 + S ),, fe_ _ A 2 2_ 

A - 2 


(26) 


where C is constant, and under the initial conditions assumed, for example, at t = 0, 
y = 0, C becomes 




y „ 2 /i„ ' 2 

° + |— - 1 1 - 2 


2 „ 2 




s 




s 


(26a) 


The extremals are then found to be 


^extremal 


-1 ± 1 / 1 + 2 


oor c 


\ 2 _ 

- 1 c 


- 1 


2C 


(27) 


The positive value is apogee and the negative value is perigee as measured from orbital 
altitude. It is felt that this relation is simple enough for a pilot to use in a manual 
computation. 


TEST CASES 

In order to test the properties of the solutions, the exact differential equations, 
first-order solution equations, and second-order solution equations were programed on a 
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digital computer. The test orbital conditions were taken as those for ejection from a 
200-kilometer circular orbit about the moon. The period of the orbit is 7600 seconds. 
Two particular trajectories were studied in detail. These trajectories were a synchro- 
nous orbit and a Hohmann transfer to the vicinity of the surface. The latter trajectory 
will henceforth be referred to as the Hohmann case. Both trajectories had a perigee 
point 20 kilometers above the surface of the moon. Such orbits are suitable either for 
reconnaissance or for landing, and both orbits are of interest for the Apollo mission. 

Synchronous Case 

The chief advantage of selecting the synchronous orbit is that since it returns to 
its initial relative position after one revolution, the investigator is in a position to 
interpret errors in the trajectories in terms of the particular terms in the solutions to 
which they are due. Trajectories were run for two separate synchronous cases: one in 
which there was no out-of -plane velocity component and one in which the initial out-of - 
plane velocity component was 10 meters per second. The latter corresponds to an orbit 
plane change of 0.36°. This out-of -plane velocity was used to show that coupling in the 
z-direction is very slight. In fact, it can be seen in figures 5 and 6 that there is no 
observable difference between the two figures to the accuracy of the plots. It can also be 
seen that the second-order solutions predict the time of occurrence of apogee and perigee 
much better than do the first-order solutions. 

Ejection was in an upward direction so that apogee is reached at approximately one- 
quarter orbit and perigee at three-quarters orbit. Time histories showing the main 
results of this study are shown in figures 5 and 6. Particular points of interest are: 

(1) The error at perigee in these cases for the second-order solution is 5 kilo- 
meters. The second-order solution is therefore an improvement over the first-order 
solution where the error is 20 kilometers and seems to indicate that the equations will 
indeed be useful as a prediction method for lunar landing or reconnaissance vehicles. 

(2) The second-order y-equation is in error timewise at the end of one orbit by 
200 seconds. Since the period for this orbit is about 6800 seconds, this error is con- 
sidered to be relatively large. The x-equation is in error timewise by the same amount. 
The reason for this time error is found in the nature of the anharmonic oscillator equa- 
tion employed in solving for y. It can be shown that when the equation is restricted to 
second-order terms, the period is not a function of amplitude. If terms of higher order 
had been carried, the period would be a function of amplitude and the time error would 
be considerably smaller (ref. 4). 
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(a) x-direction. 

Figure 5.- Time history for a 200-kilometer synchronous lunar orbit showing exact, first-order, and 

second-order solutions. 



(b) y-direction. 
Figure5.- Continued. 
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Synchronous case 
x 0 = 6.882 m/sec 
y Q = 147.8 m/sec 
z 0 - 10.0 m/sec 


Second-order solution 

Exact solution 

First-order solution 



Time, t, sec 


(c) z-direction. 
Figure 5.- Concluded. 
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(a) x-direction. 

Figure 6.- Time history for a 200-kilometer synchronous lunar orbit showing exact, first-order, and 
second-order solutions with no initial out-of-plane velocity component. 



(b) y-direction. 
Figure 6.- Concluded. 


Figure 7 is a plot of the variation of x with y for the synchronous orbit of fig- 
ure 5. Hence, this shows the position of the ejected vehicle as it would be seen from the 
ejecting vehicle if the out-of-plane motion is disregarded. It is to be observed that in 
spite of the time error in both x and y, the spatial agreement between the exact and 
approximate second-order equations is very good. 



800X10 3 


Figure 7.- Variation of x with y for a 200-kilometer synchronous lunar orbit showing the respective positions of the landing vehicle as seen 

from a vehicle in reference orbit. 
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Hohmann Case 


Figures 8, 9, and 10 show essentially the same information as figures 5, 6, and 7, 
but for the Hohmann case. First-order data have not been included in these cases. The 
agreement between exact and second-order solutions is better for the y-motion than it is 
for the synchronous orbit since the ejection velocity here is much smaller. The data 
were carried for two complete orbits in order to show more clearly the nature of the 
error buildup which occurs at the end of the first orbit in the x-equation. It is observed 
that the approximate solution is better over certain portions of the trajectory than over 
others and that the error is most significant toward the end of each orbital period. This 
result has also been obtained by the authors of reference 3 for a rectangular coordinate 
system and suggests that the cause of the error is something which is retained by both 
approaches. Figure 10 demonstrates a gradual departure from the exact solutions for 
each successive orbit because of the secular term of equation (8). 

Other Launch Angles 

In figure 11, a series of trajectories is shown for different ejection angles <p 
spaced at 20° intervals. The ejection speed is the same as for a synchronous orbit. 

None of these are synchronous orbits, however, because of the direction in which ejection 
takes place. For the sake of clarity, only the exact and second-order solutions are 
shown. 

Figure 12 contains the same information as figure 11, but the speed is the same as 
that for a Hohmann transfer. Because of the lower ejection velocity, these trajectories 
are in much better agreement than those of figure 11. 

LIMITATIONS ON APPROXIMATE SOLUTIONS 

In figure 11 no approximate solutions are presented over a range of ejection angles 
from -120° to +120°. In addition, the solutions which are presented for ±110° are con- 
siderably in error. A fundamental assumption in solving the anharmonic oscillator equa- 
tion is that the first-order term is large in comparison with the second-order term. For 
the relatively high ejection speed considered for these cases, the assumption is violated 
over this range of ejection angles. A similar situation is not encountered for the lower 
launch speed of figure 12. Mathematically, the reason for this difficulty is that for these 
trajectories, p2 becomes progressively small in' relation to the second-order term and 
then imaginary. When p2 is nearly as small as the second-order term, the solutions 

O 

are in error, and when p becomes imaginary, p is undefined and no solutions exist 
at all. 
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(b) y-direction. 
Figure 8.- Continued. 
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Figure 8.- Concluded. 
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(b) y-direction. 
Figured- Concluded. 



Figure 10.- Variation of x with y for a Hohmann transfer from 200-kilometer trajectory showing successive positions of the vehicle as seen by an 

observer in the orbiting vehicle in the reference orbit. 
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CONCLUDING REMARKS 


Second-order solutions to the equations of relative motion for two bodies in orbit 
about a planet have been presented in this paper. Also, since both a Hohmann transfer 
and a synchronous transfer from orbit to the letdown phase of a lunar landing maneuver 
are considered important with regard to current projects to explore the moon, these two 
trajectories have been studied in detail. Various other trajectories having the same 
initial speeds as these two cases but with different initial directions of motion have also 
been investigated. As a result of these studies, the second-order solutions have been 
found to be more accurate than the corresponding first-order solutions for the same types 
of motion. It has been established that errors are small for the Hohmann transfer, or 
indeed for any kind of trajectory at Hohmann speed. However, for a synchronous transfer 
the errors are somewhat larger. A certain range of trajectories at synchronous speed, 
but in nonsynchronous directions, are found to give large errors which are due to viola- 
tion of the conditions under which the second-order solutions are derived, and within this 
range is a second range in which no solutions are possible at all. These conditions are 
not thought to arise in any practical situation, however. It is therefore concluded that 
the second-order solutions may be used as a part of a guidance system where prediction 
of the future position of the space vehicle is required. 

In addition, a simple technique for computing apogee and perigee altitudes is illus- 
trated for the landing vehicle transfer trajectory when the initial separation conditions 
are known. Because of its simplicity, this technique may be used in manual guidance as 
well as automatic guidance. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., March 29, 1966. 
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DERIVATION OF THE EQUATIONS OF MOTION 


Differential Equations 

The coordinate system in this development is shown in figure 1. The Lagrangian 
is set up formally in a cylindrical coordinate system centered on the planet and then con- 
verted to shell coordinates before taking the appropriate derivatives. The Lagrangian 
in cylindrical coordinates is given by 

L = i m fr 2 + r 2 (e - uj) 2 + z 2 l + m I* (At) 

2 |/ r 2 + z 2 

This equation is converted to shell coordinates by means of the following substitutions: 


Hence, 


y + r_ = r 


x= r g 0 


z = z 


(A2) 


• • A 

y = r 

x = r g 0 > (A3) 


Then the Lagrangian becomes 


L = 




+ z* 


+ 




(A4) 


It can be seen that equation (A4) is cyclic in the x-coordinate ; thus a first integral 
of the equation of motion in the x-direction is found immediately. If the appropriate 
derivatives of the Lagrangian are taken and use is made of the exact orbital expression 
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the equations of motion become 




y - (y + r s) 


l(?r“) '“ 2 [( 1+ ^) 2+ (^t 


n -V2] 


= 0 


(A5) 


(A6) 


(A7) 


2 

Z + O) z 


For convenience, let 


K) 2 + (i ) 2 


>1 -3/2 


= 0 


X = — 
r s 


Y= -Z- 
r s 


Z = — 


T = Ci)t 




(A8) 


(A9) 


The resulting equations are 


(X - 1) = 


K 


(1 + Y) 2 

Y - (Y + 1) {(X - l) 2 - |jl + Y) 2 + Z 2 J } = 0 

r -, -3 / 2 

Z + Z [(1 + Y) 2 + Z j =0 


(A10) 


(All) 


(A12) 


where the dot over the symbols now refers to derivatives with respect to t. The 
approach taken in obtaining the solutions is the following: Since equation (All) contains 
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* 2 2 
terms in X but not in X, if Z is assumed to be small in relation to (1 + Y) and 

can hence be neglected, it is possible to cast equation (All) as an equation in Y and Y 

by direct substitution of equation (A10). This equation can then be solved approximately 

for Y as an explicit function of time. This solution may then be used in the solution of 

equations (A10) and (A12). 

In order to establish a value for K, assume that at 


T = 0 (t = 0) 


x = x o 

Y = 0 


(A13) 


Then K = ^X Q - 1^. Substituting equation (A10) into equation (All) and neglecting the 
out-of -plane term Z results in 


Y - 




(1 + Y) 3 (1 + Y) 2 


(A14) 


If the forms of 


1 

(1 + Y) n 

where n = 2 or 3 are expanded, and the terms of the second order and lower are 
retained, equations (A10), (All), and (A12) become 

X + 2KY - (K + 1) = 3KY 2 (A15) 

Y + (c* 2 )y - (|3 2 ) = (-X)Y 2 (A16) 

Z + Z = 3YZ (A17) 

where 

MV 1 ) 

a 2 = 3K 2 - 2 
j3 2 = K 2 - 1 
X = -6K 2 + 3 


(A 18) 
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First-Order Solutions 

The first-order solutions to equations (A15), (A16), and (A17) are obtained by 
dropping the second-order terms (setting the terms on the right-hand side equal to zero). 
The Y-equation is solved by inspection. Its solution is 

fl 2 

Y = a* cos (cut + e ') + £_ (A19) 

or 

where primes are used to distinguish first-order equation integration constants from 
those of the second order. 

The X-equation becomes, upon replacing K + 1 with its equivalent by defini- 
tion X Q 


X = X Q - 2KY 


(A20) 


Then substituting equation (A19) into equation (A20) and integrating, the X-equation is 


X = 



^ Ka ' sin(o!T + e ') + 
a 


2Ka' 

a 


sin e ' 


A solution to the Z -equation is 


(A21) 


Z = Z Q sin t 


(A22) 


where Z = 0 at r = 0 and Z = Z Q . 

To compute the integration constants from the initial conditions, set Y = 0 at 
t= 0 and Y = Y Q . Then solve for a* and e'. This solution can be found in a 
straightforward manner; but, for purposes of making the second-order equivalent devel- 
opment more understandable, it is advisable to make the substitution 


u' = a' cos e ’ 


and hence 


u' 


*s- 


(A23) 


Then 
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/* \ 2 

a' = \/u’ 2 + (— ) 

\ a j 


(A24) 


and 


e ' = cos -1 

a' 


(A25) 


The sign selection on equations (A24) and (A25) can be done empirically. For 
instance, a* can be chosen always positive and then the quadrant for e ' depends upon 
the launch angle <p. Equation (A19) can be expressed in a deterministic form without 
the epoch angle as 


B 2 Y o • - 2 

Y = - cos ar + — sin ar + 

a 2 « 


(A26) 


Expanding cos(q!t + e ') in equation (A19) yields 


Y = (a' cos e')cos ar - (a* sin e ')sin a.T + &=■ 

or 


(A27) 


Comparison of equations (A26) and (A27) gives 


> (A28) 


(A29) 

Since a is a frequency which is a physical quantity, a negative a has no 
meaning. Hence, with a positive, the quadrant for e ’ is determined by the respective 

O • 

signs in the numerator and denominator. Also B is strictly a function of X Q and 
always has the opposite sign (/3 2 = X Q (x o - 2), since X Q « 2^. It then follows, once the 
quadrant of e ' is known, that equations (A28) give the sign to be selected for a' in 
equation (A23). It is found that a' is always positive. These characteristics for all 
possible ejection quadrants are given in table I. Table I relates the e* -quadrant and the 
^-quadrant for the first-order solutions. 


Hence, 


a' cos € ' = 


-5 


a T sm e - - 


a 


tan e T = 


. -*o/ a 
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Second -Order Solutions 

Solution for Y . - Equation (A16) can be expressed as 

Y + a 2 Y + AY 2 = /3 2 (A30) 

It is advantageous to transform equation (A30) into the form (to eliminate the constant 
term) 

V + p 2 V + AV 2 = 0 (A31) 

This transformation can be accomplished by a change in variable Y = V + M which 
results in 

V + (a 2 + 2AM )v + AV 2 = /3 2 - AM 2 - a 2 M (A32) 

Comparing equations (A30) and (A31) shows that 

;3 2 - AM 2 - a 2 M = 0 (A33) 

and 

a 2 + 2AM = p 2 (A34) 


which can also be written as 


p 2 = ^-15K 4 - 24K 2 - 8 


(A3 5) 


Solving equation (A33) for M yields 


M 




1/2 


2A \ A 


4A" 


(A3 6) 


which can also be expressed in terms of p^ as 

2 2 
M = - — — ^-2— 
2A 


(A37) 


Equation (A31) may be recognized as the equation of a one -dimensional anharmonic 
oscillator. An exact second integral of this equation may be obtained if desired; however, 
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as this solution is an elliptic integral, the results are not very useful for solving the 
X- and Z -equations. For this reason an approximate perturbation solution is sought. 

A more complete treatment than that given here using higher order terms can be found 
in reference 4. The solution V is obtained as the sum of a first-order solution 
and a second -order correction term V 2 . Let 

V = V x + V 2 (A38) 

The first-order solution to this equation (by assuming p^ » XV) is 

Vj = a cos (pr + e) (A39) 

from which equation (A38) becomes 

V = a cos(pT + e) + V 2 (A40) 

Applying these relations to equation (A31) gives 

-p a cos (pr + e) + V 2 + p a cos (pr + e) + p V 2 + Xa cos (pr + e) 

+ 2aV 2 cos (pr + e) + Vg 2 = 0 (A41) 

After some cancellation and rearrangement of terms, 

V 2 + p 2 V 2 = -Xa 2 cos 2 (pr + e) - 2XaV 2 cos(pT + e) - XVg 2 (A42) 

Omitting terms of higher order than the second (last two terms on right), results in, 

V 2 + p 2 V 2 = -Xa 2 cos 2 (pT + e) 

= - — Xa 2 - — Xa 2 cos 2 (pr + e) (A43) 

2 2 

Solving the inhomogeneous linear equation in the usual way yields 

2 2 

V 9 = - + -—g - COS 2(pt + e ) (A44) 

Z 2p 2 6p 2 
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Hence, the second-order solution for V is 

\ a 2 j ,„2 

V = a cos (pr+e) =•+ — 5- cos 2(pr + e) (A45) 

2 p z 6 p z 

The constants a and e are integration constants which depend on initial condi- 
tions. This solution is, of course, limited to cases where the first-order assumption is 
approximately valid; that is, 

XV « p 2 


or 


A(Y - M) « p 2 


It can be seen from numerical solution that this relationship is usually the case since Y 
seldom exceeds 0.1 for cases which are physically practical and M will be small pro- 
vided the departure velocity is small enough. 

Converting equation (A45) to Y notation yields, 

Y = a cos(pT + e) - jl - | cos 2 (pr + e)J - — ^ (A46) 


The last term on the right is part of a small correction term which takes account 
of the change in energy and hence the change in average altitude between the original 
circular orbit and the new orbit into which the vehicle is launched. 


Solution for X .- The second-order solution for X is obtained by substituting 
equation (A46) for Y into equation (A15) for X and then integrating. The differential 
equation is 


X = X Q - 2K a 


2 2 2 2~1 P 

a cos (pr+ e) + ^ cos 2(pt + e) - ^ - - £L+ + 3K a 2 cos 2 (pr + e) 

6 2X 2AJ 



3 2 2 2 4 

^^-cos(pr+ e) - — ^cos(pr+ e) + cos (pr + e) - cos 2(pr + e) 

3 p' 2 x x 6 p 4 


2 2 2 

— » cos 2 (pr + e ) + cos 2 (pr + e 

6p z 6 


) 


(A47) 
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¥ 


Integrating equation (A47) with respect to r with A as an arbitrary constant of inte- 
gration yields 

X = A + Bt + C sin(pT + e) + D sin 2 (pr + e) + E sin 3 (pr + e) + F sin 4 (pr + e) (A48) 

where 


A = -(C sin e + D sin 2e + E sin 3e + F sin 4e) 

\2 


B = K 


K Ac?- P 2 ) , 3(o^ - p2 f , <A , 3o?a 2 , 19 X 2 a 4 

~ ° p2 2p2 “ ' 


K 


4X" 


24 p 4_j 


C = -K 


2a x 3a(o? - p 2 ) | 5Xa 3 

P Ap 2p 3 J 


\P 4p3 6p3 4p3 , 


E = 


F = 


KAa J 

6p 3 

KA 2 a 4 

96p 5 


(A49) 


Terms E and F have been found to be negligible in all cases tested, but are included 
here for the sake of completeness. The constant A is evaluated on the usual assump- 
tion that at t = 0, X = 0. 

Solution for Z . - The solution for Z is obtained by replacing the right-hand side 
of equation (A17) by the first-order term for Z and the appropriate second-order term 
for Y to obtain a time -dependent right-hand side. In doing so, of course, terms are 
carried which are higher than the second order of smallness. However, if only second- 
order terms are retained, new end conditions corresponding to a and e would have to 
be computed as a peripheral calculation; this calculation is thought to be unnecessary. 
Therefore, equation (A17) becomes 


Z + Z = 3 (z Q 




(A 50) 
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Let 


Then 





(A51) 


Z + Z = y cos (pr + e)sin r + 77 cos 2 (pr + e)sin r + k sin r 


(A52) 


The solution of this equation is 


where 


Z = Z Q sin r + ijyi// cos (pr + e) + rjA cos 2(pr + e) + k] sin r 
- |y i// sin(pT + e ) + ^ sin 2(pr + e) + pcj cos t 


A = 1 + p^ + p^ 

9 4 

* = i + el+ei 

4 16 

c = - — sin e - 2A sin 2e 
p 2p 




(A53) 


(A54) 


Evaluation of the Integration Constants a and e 

In order to make effective use of equation (A46), it is necessary to compute values 
for the two integration constants a and e . It is found from experience that when 
synchronous or very nearly snychronous orbits are under consideration, one can obtain 
excellent results by using the values which are obtained from first-order theory. For 
even moderate departures from the synchronous condition, for instance, on the order of 5° 
in launch direction, however, this is not the case, especially for e . Therefore, it will be 
necessary in most cases to evaluate the integration constants. A method for evaluating 
these constants is outlined in the following. 
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Differentiating equation (A46) and applying the initial conditions yields the simul- 
taneous set of transcendental equations 


0 


a cos e 


' i? ( l ■ 5 cos 2e ) 


-a? 


2 X 


• 'i o ^ 

Y = -pa sin e - sin 2e 
o K 3p 


(A55) 

(A56) 


where a and e are the two unknowns. These two equations can be solved by elimi- 
nating e in favor of the two variables u and a where 


u = a cos e 


After substitution and some manipulation, 


and 



2 Aa 

3 p 2 



(A57) 



Solving this set by eliminating 


gives a quartic equation in u 





3 P 2 (« 2 - P 2 ) , 

4 A 2 



0 


(A58) 


(A59) 


• • 

In equation (A59) it can be seen that u is a function of both Y Q and X Q (through 
a , p, and A) instead of just X Q as was the case under first-order theory (eq. (A23)). 

As a result, the special case where u = 0 will no longer occur independent of Y . It 
is seen that when u = 0, the following relation is obtained: 


3 PV-P 2 ) 

4 A* 



(A60) 


or 



(A61) 
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This relationship is shown in sketch 1. It 
can be seen that under first-order theory, 
this curve degenerates into a straight verti- 
cal line. 

This curve forms the boundary for the 
selection of the quadrant for e since the 
sign on u along with the positive selection 
for a determines the sign of cos e. 

Since the second-order term in equa- 
tion (A46) is a function of cos %(pr + e), the 
quadrant is not ambiguous and all four 
quadrants must be used. Table IE summa- 
rizes the second-order sign selection. In 
terms of ejection angle <p, equation (A61) 
becomes 



tancM— 


X, 


4 A 2 


X, 


(A62) 


Equation (A62) relates <p and X 
specifies the combinations of 

Figure 3 is a plot of x Q 


for the condition that u = 0 or in other words, it 
<p and X Q at which e changes quadrants. 

as a function of <fi for u = 0. The figure can be used 


therefore to determine for the known value of x Q and <p, the proper quadrant for e. 


LIMITATIONS ON THE APPROXIMATE SOLUTIONS 


It will be recalled that one of the conditions implicit in solving the equation of an 
anharmonic oscillator was that p^ » XV where V was defined as V = Y - M. Hence, 
in order for the solution to be valid, 


p 2 » -X(Y - M) (A63) 

where Y is inherently nondimensionalized in such a way that for all practical purposes 
the term Y in the expression can be neglected. Thus, the criterion which sets a limit 
on the validity of the solutions in actual practice is 

p 2 » -XM (A64) 
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Figure 13.- Comparison of the magnitudes of p2 and -AM for Hohmann and synchronous speeds. 

For synchronous speed at ejection angles between 110° and 250°, equation (A64) is not 
satisfied as can be seen in figure 13. Hence, one would not reasonably expect equa- 
tion (A46) to describe the physical situation. That this is actually the case can be seen 
by looking at figure 11 for ejection angles of 110° and -110° (250°). It is observed that 
for these cases the error is considerable. At ejection angles of slightly larger magni- 
tude than this value, p 2 becomes imaginary (p 2 = \[ 4A/3 2 + a^). Hence, a value for p 

is undefined. Inspection of equation (A46) shows that p is the frequency term. There- 

2 

fore, no solutions would be expected to exist at all with p imaginary. This situation 

occurs for those cases in figure 11 for which no solutions are shown. For the Hohmann 

o 

case, equation (A46) is satisfied for all ejection angles and p remains real; hence, 
breakdown does not occur. 
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CONVERSION OF DATA FROM RECTANGULAR TO SHELL COORDINATES 


If initial measurements are taken in an x',y',z' rectangular coordinate system, 
the conversion to an x,y,z shell coordinate system can be derived as follows where all 
terms are dimensional. 


x’ 



From sketch 2, where it is assumed, as before, 
that r g is measured toward the center of the planet, 


tan 8 = 


r s + y' 


but 8 = — . Therefore, 
r 

A s 


x = r tan -1 — — — 
s r s + y f 


and from the large right triangle in the sketch 
y = + /x' 2 + (y' + r g ) 2 - r g 


(Bl) 


(B2) 


where the positive sign is selected for the radical because r and the radical are of the 

s 

same order of magnitude. 

Equations (Bl) and (B2) serve to express spatial positions x and y in shell 
coordinates in terms of x' and y' in rectangular coordinates. Of course, the out-of- 
plane measurement is trivial since 


z = z' 


(B3) 


The velocity terms are obtained by differentiation of equations (Bl), (B2), and (B3) 
with respect to time to get 


x = r. 


(r g + y’)x’ - x'y'" 


x' 2 + 


( r s + y') 


t\2 


. x’x' + (y + r )y’ 

y = r . = N =t — 

i /* ,2 + ( y + r s ) 2 


Z = Z 


(B4) 

(B5) 

(B6) 
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TABLE I.- FIRST-ORDER QUADRANT AND SIGN SELECTION 



a' 

+ 

+ 

+ 

+ 


(j) quadrant 

1 

2 

3 

4 


TABLE E.- SECOND-ORDER QUADRANT AND SIGN SELECTION 


Ejection quadrant 

Y 0 

u 

Sign 
on a 

Quadrant for e 

0 = 0 = 17 

Y„>0 

u > 0 

+ 

e in fourth quadrant 

VII 

■©- 

VII 

o 

o 

A 

O 

•>* 

u < 0 

+ 

e in third quadrant 

i7 = 4> = 2ir 

Y o <0 

u < 0 

+ 

€ in second quadrant 

17 = <p = 2ir 

- 

Y < 0 

u > 0 

+ 

e in first quadrant 
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